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Abstract 

We analyze the spatially semidiscrete piecewise linear finite element method for 
a nonlocal parabolic equation resulting from thermistor problem. Our approach 
is based on the properties of the elliptic projection defined by the bilinear form 
associated with the variational formulation of the finite element method. We assume 
minimal regularity of the exact solution that yields optimal order error estimate. 
The full discrete backward Euler method and the Crank-Nicolson-Galerkin scheme 
are also considered. Finally, a simple algorithm for solving the fully discrete problem 
is proposed. 
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1 Introduction 



We study the numerical approximation by the finite element scheme of the 
nonlinear problem 

^ - V.{k{u)Vu) = X- ^^""^ ^„ in Qx]0; t[ 

( f(u) dx 



u = on dilx]0; t[, 
u{0) = Uq in Q, 
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where is a bounded domain in M^, t is a positive fixed real, / and k are func- 
tions from R to M satisfying tlie liypotlieses (HI) — {H2) below, A is a positive 
parameter and V denotes the gradient with respect to the a;— variables. The 
time evolution model (1) describes the temperature profile of a thermistor 
device with electrical resistivity /, see [3,4,5,8,9,10,14,15,17,19]; the dimen- 
sionless parameter A can be identified with the square of the applied potential 
difference V at the ends of the conductor. The system (1) has been the subject 
of a variety of investigations in the last decade. Existence of weak solutions 
to problems related with the thermistor problem is proved in [16], where the 
mathematical treatment of this system apparently appears for the first time. 
In [11] the problem (1) for the special case /c = 1 is considered, and then a 
backward Eulcr time-semidiscretization method for the approximation of its 
solution is proposed and analyzed. In this paper we propose a finite element 
method to construct numerical approximations of the solutions of problem (1) 
for the case when k is different from the identity. The formulation of the finite 
element method is standard and it is based on a variational formulation of the 
continuous problem. There is a vast literature on finite element methods for 
nonlinear elliptic and parabolic problems. For example, we mention the work 
[7] on the porous media equations, which are similar to the Joule heating prob- 
lem [1] . Compared to a standard semilinear equation, the main challenge here 
is the nonstandard nonlocal nonlinearity on the right-hand side of the partial 
differential equation (1). 

On the other hand, error bounds are normally expressed in terms of norms of 
the exact solution of the problem. It is well known that the required regularity 
of the exact solution can be attained by assuming enough regularity of data, 
sometimes supplemented with compatibility conditions, see [2,13,20]. We then 
use sufficient conditions in terms of the data of the problem and its solution 
u that yield error estimates (see hypotheses (H1)-(H3) below). 



2 Main results and organization of the paper 

We denote by (•, •) and || • || respectively the inner product and the norm in 
— -L^(Q), by \\-\\s the norm in the Sobolev space H^{Q), by c some generic 
positive constant which may depend upon the data and whose value may vary 
from step to step. In Section 3 we study spatially semidiscrete approximations 
of (1) by the finite element method. The approximate solution is sought in 
the piecewise linear finite element space 

Sh = Sh{n) = {x e C(Q) : x/e Zmear, Ve e T^; x/sn = 0} , 

where {Th}h is a family of regular triangulations of ^2, with h denoting the 
maximum diameter of the triangles of T^. As a model for our analysis we 
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first consider the corresponding semidiscrete Galerkin finite element method, 
which consists in finding Uh{t) £ Sh such that 

{uh,t,X) + Vm/j, Vx) = 7 r2ifiuh),x), 

(^/K)dx) (2) 

Uh{0) = uoh , 

Vx € Sh, t E J = (O.t), and where Uoh E Sh is a given approximation of 
uq. Similar discretization techniques have been analyzed for various linear and 
nonlinear evolution problems (cf. e.g. [13]). This method (2) may be written as 
a system of ordinary differential equations. In fact, let be the standard 

pyramid basis of Sh- Write Uh{x,t) — Y,f=iC(j{t)4>j{x), where {oij)i<j<Nh are 
the real coefficients to be determined. Then, (2) can be written as 

Aa'{t) + B{a)a{t) = f{a), teJ, a(0) = 7, 

where 7 is the vector of nodal values of uoh, /(o;) — {fi{(^), ■ ■ ■ , fN^i'^))'^ with 



A 



and A = {ajk)i<j,k<Nh and B{a) = {hjk{(y))\<j,k<Nh are, respectively, the cor- 
responding mass and stiffness matrices: 



We shall assume the following general assumptions on the given data: 

(HI) / :]R^]Risa locally Lipschitzian function and f{u)>a>Q for all 

M e R; 

(H2) A; is a twice derivable function verifying: there exist positive constants 

/ci, k2 and c such that < /ci < k{u) < k2, \k'{u), k"{u)\ < c; 
(H3) u e L°°{0,T,H^{Q)nW^'°°{n)) and uq G H^{n). 

It is shown in [6, Theorem 2.1] that the regularity assumption {H3) is satisfied 
if the data is smooth and compatible. The matrix A is always definite positive. 
Further, hypotheses (HI) and (H2) assure that the matrix -8(0;) is also positive 
definite. Assumption (H3) is used in order to prove Lemma 1 and then to show 
a o{h'^) error estimate. It is useful to introduce the interpolation operator 
Ih : C{fl) Sh defined by 

IhV = J2v{Pj)(l)j{x), 
i=i 
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where {Pj}jjli are the interior vertices of T^. For completeness, we assume the 
following standard interpolation error estimates : for v G ni7o(0) there 

exists a positive constant c > such that 

Whv — v\\ < c/t^ll't'lh and \\V{IhV - v)\\ < ch\\v\\2, 
holds. We also assume the property [18]: for some integer r > 2 and small h, 

inf {11^; - xll + h\\Viv - x)||} < ch'\\v\\s, for 1 < s < r, (3) 

when V G H'^CiHq. We finally suppose that the family of triangulations is 
such that the inverse estimate [18] 

||Vx||<c/i-^||x|| Vxe5, (4) 

is satisfied. In the existing literature the error estimates for the finite element 
method are normally expressed in terms of norms of the exact solution of the 
problem and are usually derived for solutions that are sufficiently smooth (cf. 
e.g. [13]). To estimate the error in the semidiscrete problem (2) we split the 
error as Uh — u = {uh — Uh) + {uh — u) = 9 + p, where Uh denotes the standard 
elliptic projection in Sh of the exact solution u defined by 

(k(u(t))V(uh-u),Vx)^0, ^x^Sh, t>0. (5) 

It is well known (see [12,18]) that if the regularity hypothesis (H3) for u holds, 
then Uh has the following approximation properties: 

Lemma 1 Under the regularity hypotheses (H1)-(H3), one has 

\\p{t)\\+h\\vpm<c{u)h', 
\\ptm + h\\vptm<c{u)h', 

where c{u) is a constant independent of t & J. 

Lemma 2 LetUh be defined by (5). Then, ||Vitfe(t)||Loo < c{u), t E J. 

For the finite clement method (2), if Uoh is chosen such that \\uoh — ^o|| ^ 
c/i^||Mo||2, we prove in Section 3 an error estimate of the form \\uh{t) —u{t)\\ < 
c{u)h^ (Theorem 3). The corresponding estimate for the gradient is also proved 
(Theorem 4). In Section 4 we show that our approach for the semidiscrete 
Galcrkin finite element method also applies to fully discrete schemes. We con- 
sider the backward Euler method for the discretization in time of (1): letting 
r to be the time step, the approximation in Sh of u{t) at t = t„ = nr, 
i.e. [/" = J2f=iOi]4>j, where {c(^)i<j<Nh are the unknown real coefficients. 
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ljn_un-l 



, n — 0,1,2, ... , the numerical method is defined by 



T 



{dnU^,x) + {k{unvu\Vx) 



A 



{j^fiUn)dx) 



(6) 



For this scheme we prove (Theorem 5), under the same regularity requirements 
as in Section 3, that 



In order to obtain higher accuracy in time, in Section 5 we investigate an 
alternative way to obtain an o{h^ + r^) error bound using the basic Crank- 
Nicolson-Galerkin scheme: applying the techniques of sections 3 and 4 we prove 
(Theorem 7) an error estimate of the form 



Finally, in Section 6 we propose a simple algorithm for solving the fully discrete 
problem. 



3 Semidiscrete problem 

In this section we obtain an error estimate for u and the associated estimate 
for the gradient. The proofs use a sphtting of the error based on the eUiptic 
projection Uh (5). We may define the semidiscrete problem on a finite interval 
J — {0,t\ of time. 

Theorem 3 Let u andu^ he the solutions of (1) and (2), respectively. Then, 
under the hypotheses (H1)-(H3), we have: \\uhit) — u{t)\\ < c{u)h'^, for t E J, 
provided that \\uoh ~ Uo\\ < ch?. 



PROOF. Owing to Lemma 1 and to the decomposition of the error as sum of 
two terms Uh — u — {u^ — Uh) + {uh — u)=9 + p,\i suffices to treat 9 — Uh — Uh- 



U"" - U{tn)\\ < C{u){h^ + T), for tn G 



J^{0,t). 



C/" - u{tn)\\ < c{u){h^ + r"), for i„ e 



J=(0,i). 
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We have from the equations satisfied by Uh and Uh that 

= {uh,uX) + ik{uh)yuhyx) - {uh,uX) - {k{uh)^Uh,Vx) 

^ {J f{uh)Y ^^^^^^'^^ ~ ~ ~ (^NVwfe, Vx) + {{k{u) - k{uh))Vuh, Vx) 

^ r(/K), X) - {pt, X) - {ut, X) - {k{u)\/u, Vx) + {{k{u) - k{uh))VUh, Vx) 



(//K)) 



2 ' 



^ r(/K), X) - 77^^(/H' X) + {{kH - k{uh))VUh, Vx) - (pt, x) 



^ (/K) - / W), X) + {jT^-Tyi - TT^^ im,X)) 



+ {{k{u) - k{uh))Vuh, Vx) - {pt, X) 

^ 7777 — v^i-fi^h) - /(«)), X) + {{k{u) - k{uh))Vuh, Vx) - (pt, x) 

+ (J/K)W/M)- (X(/M-/M).ix) (/^(/K) + (/(„), X). 

Thus, setting x = using the hypotheses (H1)-(H3), Lemma 1, and Young's 
inequahty, 

^^ll^ll' + k,\\ver < c{\\u, - u\m\ + liv^^ii) + iip.ii ii^ii) 

<k,\\ver+c{\\er+\\pr+\\ptr)- 

By integration, we get \\e{t)f < ||^(0)p + c/o(||^p + \\pf + \\ptf)ds and it 
follows by Gronwall's Lemma that 

\\e{t)r<c\mr+c[{\\pr+\\ptr)ds 

Jo 



or 



||6'(0)|| < \\uoh - UqW + \\uh{0) - Uo\\ < \\uoh - UoW + ch'^\\uo\\2- (7) 
We then get the desired conclusion: 

\\0{t)\\ < c\\uoh - UoW + c{u)h^ < c{u)h^. 



We now derive for the standard Galerkin method, from the weak formulation 
of the parabolic problem and using the inverse property, the optimal order 
error estimate for the gradient. 

Theorem 4 Letu anduh he, respectively, the solutions of (1) and (2). Under 
hypotheses (H1)-(H3), if Ugh is chosen such that \\uoh — Uo\\ < ch'^\\uo\\2, then 

\\Vuh{t) - Vu(t)\\ < c(u)h, for t e J. 
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PROOF. We have 

\\Wuh{t) - Vu{t)\\ < ||VK(t) - + ||Vx - Vu{t)\\ 

<ch-'\\un{t)-x\\ + \\^X-^uit)\\ (8) 
< ch-'\\uh{t) - u{t)\\ + ch-\\\x - u{t)\\ + h\\Vx - Vu{t)\\). 

By the approximation assumption (3) we know that, with suitable % e Sh, 

\\X - u{t)\\ + h\\Wx - ^<t)\\ < ce\\u{t)\\2. 

Then, using (8), we get: 

\\Vuh{t) - Vu{t)\\ < ch-'\\uhit) - u{t)\\ + ch\\u{t)\\2. 

Theorem 3 yields the intended conclusion: ||Vii/i(t) — Vii(t)|| < c{u)h. 

4 The completely discrete case 

We now turn our attention to the fully discrete scheme based on the backward 
Euler method: find f/" G Sh such that (6) holds. Existence result of (6) is a 
simple consequence of Brouwer's fixed point theorem. Here we obtain an error 
estimate for the scheme. For h, r and ^ small enough, we prove uniqueness. 

Theorem 5 Let u and U"- be solutions of (1) and (6) respectively, with UQh 
chosen such that \\uq — UQh\\ < ch"^. Under the required regularity (H1)-(H3), 
there exists a constant c such that, for tn & J and small t, we have 

W-U{tn)\\<c{u){h''+T). 

Moreover, for sufficiently small r, h and j^, there is a unique solution C/" of 
the complete discrete scheme (6). 

PROOF. We use the partitioning of the error 

- + (^/- - =9^ + Pn, (9) 

with = u{tn), U"' = Uh{tn), whcrc Uh is the elliptic projection of -u" defined 
by (5). By virtue of Lemma 1, it suffices to bound 6'„. We have for x ^ Sh that 

(9„^n,x) + WC/")V^„,Vx) 

= ac/",x) + (A;(nVC/^ Vx) - (dnU^x) - {k{U^)VlJ\Vx) 

= 7 r2(/(f^"),x) - - (5n[/" -<,X) 

{f^fmdx) 

- {k{u^)VU\ Vx) - ((^(n - k{u^))VU\ Vx) 
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and, in view of the elliptic projection (5) and the equation of the continuous 
problem (1), we can write: 



(Mn,x) + (Mt^")ve„,Vx) 

- 2(/(n,X) 



A 



r(/K),X) 



{j^f{u-)dx) {j^fiu-)dxy 

X 



{j^f{U-)dx) 



r(/(n-/K),X) 



A 



r)(/K),X) 



- {dnPn,X) - -<,X) - - k{u^'))VU\Vx). 

Choosing x — and using the fact that VC/", vJ^, C/" are bounded, we get: 

\dn\\en\? + h\\ve,,f 

< c||c/" - i^lKli^nii + ||v^„||) + + WdnU'^ - <||)||^„||. 

Hence, by Young's inequality, 

^nll^nll' + hwvej' < c{\\ej' + yj' + Wd^PuW' + l|5n«" - (10) 

Introducing the notation i?„ = ||p„|p + ||9„p„p + ~ we write (10) 

in the form (?n||^n|P + ^i||V6'„|p < c(||6'„|p + i?„), and it results that 

il-Cr)\\enf<\\en-lf + CrRn. 

For r < i, we have < iz^||^n-i||^ + f^- Using the fact that « 

1 + CT for T sufficiently small, it follows that ||6'„|p < (1 + cr)||6'„_i|p + cri?„. 
By induction, we get: 



ll^nll' < (1 + CT)"||eo||' + cr 5](1 + crr-^Rj 

n 

< c\\9of + CT ^ Rj, for tn e J. 



(11) 



We now recall that from Lemma 1 

\\p,\\ < c{u)h\ \\d^p,\\ 

On the other hand, we have 



1 /•* 



T Jt 



Ptds 



■3-1 





2 




dnU^ - ul 




-/ 




r Jt, 



< c{u)h^ 



< c{u)t. 
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It yields that Rj < c{u){h? + r)^. Taking the above estimates together with 
(11) and (7) we prove the intended bound for 0^. 

\\0n\\ < c\\uoh - uo\\ + c{u){h^ + r) < c{u){h'^ + r). 



It remains to prove the second part of the theorem (uniqueness) . Let U"' = X 
and [/" = y be two solutions of the fully discrete problem: 

{X-Y,x)+T {k{X)VX - k{Y)VY, Vx) 

At , , / A- A- 



Taking x — X — Y, we have 

\\X - Yf + T{k{X)V{X - Y), V{X - Y)) = ^J7^(/ W - fiY\X - Y) 

^ [iiww ~ uww) ^ " " ^ ^^^^^^ " " 

Thus, 

\\X - Yf + lkir\\V{X - Y)f < c\\X - Yfir + T||Vy||ioc). 

According to Lemma 2 we have 

||Vy||Loo < ||VSft||Loc + ||V6'„||loo < c + c/i~^||V6'„||. 
Taking into account the estimate for we get 

A;i||V^„||^ < C||^„_i||^ + t||^„||2 + TRn < c{h^ + T)^ 

and we deduce that 



rllVyllioc < c(t + /i' + 



2^ 

2 



Then, for sufficient small r, h and ^, we get the uniqueness of the solution of 
the complete discrete scheme. 
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5 The Crank-Nicolson-Galerkin scheme 



This section is devoted to the study of the following Crank-Nicolson-Galerkin 
scheme: 



[!nf{U )dx) 



(12) 



o{h? -\- T^) error bound - we need an auxiliary estimate. 



with U — -'"^ — . Before proceeding to the main result of the section - an 



Lemma 6 Besides the hypotheses (H1)-(H3), let us further assume 
(H4) UueL°-{0,T,H\n)). 

Then, \\Vuh^tt\\ < c(m), where Uh is the elliptic projection defined by (5). 



PROOF. Differentiating (5) twice in time we find that 

{k{u)Vuh,u,^x) = {k{u)Vuu,Vx) - 2{k'{u)Vpt,Vx) - {k"{u)Vp,Vx). 
Taking x = Uh,tt, it follows that 

ki\\^Uh,tt\\^ < c{u){\\Vutt\\ + llVptll + ||Vp||)||Vfift,«|| 
and, recalling Lemma 1 and (H4), we obtain the intended conclusion. 

Theorem 7 Letu and be respectively the solutions of (1) and (12). Then, 
under the hypotheses (H1)-(H4), we have \\U"' — u{tn)\\ < c(M)(/i^+r^), i„ e J. 
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PROOF. Splitting the error by means of the eUiptic projection as in (9), one 
gets 



(a„r,x) + (A:([/")Vr,Vx 

= (5nt/",x) + (^(r)vir, vx) - (a„t/^x) - [k{U'')VU'\Vx 



A 



{hf{Tr)dxy 

- (k VC/"-5, Vx) - (fc(ir)V^ - k VC/"-^ Vx) 



A 



(/K-^),x) 



+ A;K-^)V (^^-C7"-5^ ,Vx). 
Setting X = ^, it follows from {dJ'',T) = and (H1)-(H3) that 

,2i(/K"^)-/(t^),vr)i 
^(/(i7"),vr) 



l9irir+Hlvrir<- V 



+ 



A 



A 



+ c(ii9„t7" - i^r^ii + - II + II - t/"-^)ii) II vr II 

< c(||a„[7" - + - + l|v(F - [/"-5)||)||v^||. 

Young's inequality gives 

dp^'f < c(ii9„c/" - ^^ri' + - ^""^ir + iiv(^ - c/"~^)ir). (is) 

Estimating each term on the right hand side of the inequality (13) separately, 
we have 

iicr-ii"-5|| < ||rii + iip"ii + i|Tz"-ii"-5|| < ||rii+c(^x)(/i2+r2), (m) 

\\dJj^-uT'^\ < \\dnP''\\ + \\dnU^-ur'^\\<c{u)(h^+r'), (15) 

and by Lemma 6 



V{U -C/"-2) 



<CT \\Vuh,tt\\ ds < c{u) 

Jt„-1 



(16) 
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Inequalities (13)-(16) together show that 

ayr ||2 < c\\T\\^ + c{u){h^ + < + r-^^ + c{u){h'' + r'f 

which gives (1 - cr)||rf < (1 + ct)!!^"-^^ + c{u)T{h^ + r'^f. Then, by 
induction, we have that for smaU r 

ll^'^ll < c\\e^\\ + c{u){h'^ + T^) < c\\uoh - UoW + c{u){h'^ + T^) for i„ e J. 

Or ||iio/i — ""oil < c/i^, and then H^^H < c{u){h^ + r^), which leads, in view of 
Lemma 1, to the intended result. 



6 Method of solution 



We divide the interval O = (—1, 1) into equal finite elements: Xq = —1 < 
xi < . . . < xn = ^- Let {xj, Xj+i) be a partition of ^2 and Xj+i — xj = h = 
the step length. By 5" we denote a basis of the usual pyramid functions: 



Vj{x) 



+ on[xj^i,Xj], 
-l^ + i'^+j) on[xj,Xj+i], 
otherwise. 



We first write (1) in variational form with k = 1. We then multiply (1) by Vj 
(for j fixed). We have, using the boundary conditions, that 



/ -rr-v-j dx + VuVv-j dx 
Jn at Jn 



{lnf{u) dx) 
By the Crank- Nicolson approach we obtain 

Jn Jn Jn (^J^f{u^)dx) 

so the approximation to the function u can be written as 

AT 



i=-l 
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where the a""^^ are unknown real coefficients to be determined. Prom (17) it 
is easy to obtain the following system of (A^ — 1) linear algebraic equations: 



h ^ 2h „ h „ Ar f ((/")(•; fix 
— — rv , -\ (T -I fy , -I — — 

Using the Dirichlet boundary conditions, we have 



and 



al+^ ^ a'l+^ - {1 + h)a:^+\ 



^n+l _ n+1 



(19) 



(20) 



where c = (jv-i)/t-i - Substituting the expressions (19) and (20) into (18), we 
obtain the system of equations 



{b-{l + h)) + 2aa'l+^ 



^-(i + /.)U^ + ^< + ^M(!^,,.o 



h „ 2/i „ h „ \t fc f(u'^)vi dx . , ^ 
^ 3^6^ (/n/(^^'^)^x) 

aQ;^t^2 + («c + 6)Qi^t\ 

(21) 

where 

h T , 2/i r 

a= — , 0= 1 — -, 

6 h' 3 2/i' 

and thus a simple algorithm for solving the fully discrete problem. 
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